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Rapid developments in geographical information systems (GIS) 
continue to generate interest in analyzing complex spatial datasets. 
One area of activity is in creating smoothed disease maps to de- 
scribe the geographic variation of disease and generate hypotheses 
for apparent differences in risk. With multiple diseases, a multivariate 
conditionally autoregressive (MCAR) model is often used to smooth 
across space while accounting for associations between the diseases. 
The MCAR, however, imposes complex covariance structures that are 
difficult to interpret and estimate. This article develops a much sim- 
pler alternative approach building upon the techniques of smoothed 
ANOVA (SAN OVA). Instead of simply shrinking effects without any 
structure, here we use SANOVA to smooth spatial random effects 
by taking advantage of the spatial structure. We extend SANOVA 
to cases in which one factor is a spatial lattice, which is smoothed 
using a CAR model, and a second factor is, for example, type of 
cancer. Datasets routinely lack enough information to identify the 
additional structure of MCAR. SANOVA offers a simpler and more 
intelligible structure than the MCAR while performing as well. We 
demonstrate our approach with simulation studies designed to com- 
pare SANOVA with different design matrices versus MCAR with dif- 
ferent priors. Subsequently a cancer-surveillance dataset, describing 
incidence of 3-cancers in Minnesota's 87 counties, is analyzed using 
both approaches, showing the competitiveness of the SANOVA ap- 
proach. 

1. Introduction. Statistical modeling and analysis of spatially referenced 
data receive considerable interest due to the increasing availability of geo- 
graphical information systems (GIS) and spatial databases. For data aggre- 
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gated over geographic regions such as counties, census tracts or ZIP codes 
(often called areal data), with individual identifiers and precise locations 
removed, inferential objectives focus on models for spatial clustering and 
variation. Such models are often used in epidemiology and public health to 
understand geographical patterns in disease incidence and morbidity. Recent 
reviews of methods for such data include Lawson et al. (1999), Elliott et al. 
(2000), Waller and Gotway (2004) and Rue and Held (2005). Traditionally 
such data have been modeled using conditionally specified probability mod- 
els that shrink or smooth spatial effects by borrowing strength from neigh- 
boring regions. Perhaps the most pervasive model is the conditionally au- 
toregressive (CAR) family pioneered by Besag (1974), which has been widely 
investigated and applied to spatial epidemiological data [Wakefield (2007) 
gives an excellent review] . Recently the CAR has been extended to multivari- 
ate responses, building on multivariate conditional autoregressive (MCAR) 
models described by Mardia (1988). Gelfand and Vonatsou (2003) and Carlin 
and Banerjee (2003) discussed their use in hierarchical models, while Kim, 
Sun and Tsutakawa (2001) presented a different "twofold CAR" model for 
counts of two diseases in each areal unit. Other extensions allowing flexible 
modeling of cross-correlations include Sain and Cressie (2002), Jin, Car- 
lin and Banerjee (2005) and Jin, Banerjee and Carlin (2007). The MCAR 
can be viewed as a conditionally specified probability model for interactions 
between space and an attribute of interest. For instance, in disease map- 
ping interest often lies in modeling geographical patterns in disease rates or 
counts of several diseases. The MCAR acknowledges dependence between 
the diseases as well as dependence across space. However, practical difficul- 
ties arise from MCAR's elaborate dependence structure: most interaction 
effects will be weakly identified by the data, so the dependence structure is 
poorly identified. In hierarchical models [e.g., Gelfand and Vonatsou (2003), 
Jin, Carlin and Banerjee (2005, 2007)], strong prior distributions may im- 
prove identifiability, but this is not uncontroversial, as inferences are sensi- 
tive to the prior and perhaps unreliable without genuine prior information. 
This article proposes a much simpler and more interpretable alternative to 
the MCAR, modeling multivariate spatial effects using smoothed analysis of 
variance (SAN OVA) as developed by Hodges, Carlin and Fan (2007), hence- 
forth HCSC. Unlike an ANOVA that is used to identify some interaction 
effects to retain and others to remove, SANOVA mostly retains effects that 
are large, mostly removes those that are small, and partially retains middling 
effects. (Loosely speaking, "large," "middling" and "small" describe the size 
of the unsmoothed effects compared to their standard errors.) To accommo- 
date rich dependence structures, MCAR introduces weakly identifiable pa- 
rameters that complicate estimation. SANOVA, on the other hand, focuses 
instead on smoothing interactions to yield more stable and reliable results. 
Our intended contribution is to show how SANOVA can solve the multiple 
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disease mapping problem while avoiding the dauntingly complex covariance 
structures imposed by MCAR and its generalizations. We demonstrate that 
SANOVA produces inference that is largely indistinguishable from MCAR, 
yet SANOVA is simpler, more explicit, easier to put priors on and easier to 
estimate. The rest of the article is as follows. Section 2 reviews SANOVA 
and MCAR, identifying SANOVA as a special case of MCAR. Section 3 is a 
"tournament" of simulation experiments comparing SANOVA with MCAR 
for normal and Poisson data, while Section 4 analyzes data describing the 
number of deaths from lung, larynx and esophagus cancer in Minnesota 
between 1990 and 2000. A summary and discussion of future research in 
Section 5 concludes the paper. Zhang, Hodges and Banerjee (2009) (Appen- 
dices) gives computational and technical details. 

2. The competitors. 

2.1. Smoothing spatial effects using SANOVA. 

2.1.1. SANOVA for balanced, single- error-term models [HCSC (2007)]. 
Consider a balanced, single-error-term analysis of variance, with M% de- 
grees of freedom for main effects and M2 degrees of freedom for interactions. 
Specify this ANOVA as a linear model: let A\ denote columns in the de- 
sign matrix for main effects, and A2 denote columns in the design matrix 
for interactions. Assume the design has c cells and n observations per cell, 
giving cn observations in total. To simplify later calculations, normalize the 
columns of A\ and A2 so = Imi an d A' 2 A2 = Im 2 - (Note: HCSC nor- 

malized columns differently, fixing A'-^A\ = cuIm 1 and A' 2 A2 = cn/jvf 2 .) Then 
write the ANOVA as 

~©1 
2 



(1) 



[A!\A 2 



+ e = Ai@ 1 + A 2 @2 +e, 



where e ~ N(0, —I) with 770 being a precision, y is cn X 1, A± is cn x 
Mi, A2 is cn x M2, 0i is M\ x 1, 2 is M2 x 1, and e is cn x 1. This 
ANOVA is smoothed by further modeling 0. HCSC emphasized smooth- 
ing interactions, although main effects can be smoothed by exactly the 
same means. Following HCSC, we add constraints (or a prior) on ©2 as 
Qm-l+j ~ N(0, 1/rjj) for j = 1, . . . , M2, written as 

"01 



(2) 



0m 2 = [0m 2 xMiI-7m 2 ] 
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N(0, diag(— )), in the manner of Lee and Nelder (1996) and 



where d 

Hodges (1998). Combining 
linear model: 
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[1) and (2), express this hierarchical model as a 
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More compactly, write 



(4) 



Y=X@+e 



where Y has dimension (cn + M2) x 1 and e's covariance T is block diagonal 
with blocks T\ = ^I C n for the data cases (rows of X corresponding to the 
observation y) and T2 = diag(l/r/i, . . . , 1/t]m 2 ) for the constraint cases (rows 
of X with error term 8). For convenience, define the matrix Xr> = [A\\A2], 
the data-case part of X. This development can be done using the mixed 
linear model (MLM) formulation traditionally written as y = Xf3 + Zu + e, 
where our (1) supplies this equation and u = ©2 ~ A(0,r2). The develop- 
ment to follow can also be done using the MLM formulation at the price 
of slightly greater complexity, so we omit it. HCSC developed SANOVA for 
exchangeable priors on groups formed from components of ©2- The next 
section develops the extension to spatial smoothing. 

2.1.2. What is CAR? Suppose a map has N regions, each with an un- 
known quantity of interest (pi, i = 1, . .. , N. A conditionally autoregressive 
(CAR) model specifies the full conditional distribution of each (pi as 



where i ~ j denotes that region j is a neighbor of region i (typically defined 
as spatially adjacent), and rrii is the number of region i's neighbors. Equation 
(5) reduces to the well-known intrinsic conditionally autoregressive (ICAR) 
model [Besag, York and Mollie (1991)] if a = 1 or an independence model if 
a = 0. The ICAR model induces "local" smoothing by borrowing strength 
from neighbors, while the independence model assumes spatial independence 
and induces "global" smoothing. The CAR prior's smoothing parameter a 
also controls the strength of spatial dependence among regions, though it 
has long been appreciated that a fairly large a may be required to induce 
large spatial correlation; see Wall (2004) for recent discussion and examples. 
It is well known [e.g., Besag (1974)] that the conditional specifications in 
(5) lead to a valid joint distribution for cf> = ((pi, . . . , (p^) 1 expressed in terms 
of the map's neighborhood structure. If Q is an N x N matrix such that 
Qu = rrii, Qij = —a whenever i ~ j and Qij = otherwise, then the intrinsic 
CAR model [Besag, York and Mollie (1991)] has density 





(6) 




with 




N-G 



if a 6 (0,1) 
if a = 1. 
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In (6) r is the spatial precision parameter, tQ is the precision matrix in 
this multivariate normal distribution and G is the number of "islands" (dis- 
connected parts) in the spatial map [Hodges, Carlin and Fan (2003)]. When 
a £ (0, 1), (6) is a proper multivariate normal distribution. When a = 1, Q is 
singular with Ql = 0; Q has rank N — G in a map with G islands, therefore, 
the exponent on r becomes (N — G)/2. In hierarchical models, the CAR 
model is usually used as a prior on spatial random effects. For instance, let 
Yi be the observed number of cases of a disease in region i, i = 1, . . . , N, 
and Ei be the expected number of cases in region i. Here the Yi are treated 
as random variables, while the Ei are treated as fixed and known, often 
simply proportional to the number of persons at risk in region i. For rare 
diseases, a Poisson approximation to a binomial sampling distribution for 
disease counts is often used, so a commonly used likelihood for mapping a 
single disease is 

(7) Yi ~ d Poisson(^e w ), * = 1,...,JV, 

where \n = x' i f3 + (j)i- The x, are explanatory, region-specific regressors with 
coefficients j3 and the parameter /ij is the log-relative risk describing depar- 
tures of observed from expected counts, that is, from Ei. The hierarchy's 
next level is specified by assigning the CAR distribution to <p and a hyper- 
prior to the spatial precision parameter r. In the hierarchical setup, the 
improper ICAR with a = 1 gives proper posterior distributions for spatial 
effects. In practice, Markov chain Monte Carlo (MCMC) algorithms are de- 
signed for estimating posteriors from such models and the appropriate num- 
ber of linear constraints on the <f> suffices to ensure sampling from proper 
posterior distributions [Banerjee, Carlin and Gelfand (2004), pages 163-164, 
give details]. 

2.1.3. How does CAR fit into SANOVA? To use CAR in SANOVA, the 
key is re-expressing the improper CAR, that is, (6) with a = 1. Let Q have 
spectral decomposition Q = VDV 1 , where V is an orthogonal matrix with 
columns containing Q's eigenvectors and D is diagonal with nonnegative 
diagonal entries. D has G zero diagonal entries, one of which corresponds 
to the eigenvector — =ljv, by convention the iVth (right-most) column in 

V. Define a new parameter = V'(f), so has dimension iV and precision 
matrix tD. Giving an iV-vector a normal prior with mean zero and preci- 
sion tD is equivalent to giving = VQ a CAR prior with precision tQ. 
consists of ©A? = -j=l' N <j> = y/N </>, the scaled average of the fa, along with 
iV— 1 contrasts in <fi, which are orthogonal to -^=1jv by construction. Thus, 
the CAR prior is informative (has positive precision) only for contrasts in 
<fi, while putting zero precision on &gm = &N = "m^N^i ^ ne overall level, 
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and on G — 1 orthogonal contrasts in the levels of the G islands. In other 
words, the CAR model can be thought of as a prior distribution on the 
contrasts rather than individual effects (hence the need for the sum-to-zero 
constraint). A related result, discussed in Besag, York and Mollie (1995), 
shows the CAR to be a member of a family of "pairwise difference" pri- 
ors. This reparameterization allows the CAR model to fit into the ANOVA 
framework, with @gm corresponding to the ANOVA's grand mean and the 
rest of 0, ®R eg , corresponding to V^'(f>, where TA _ ) is V excluding the col- 
umn -j=ljv, consisting of N — 1 orthogonal contrasts among the ./V regions 
and giving the N — 1 degrees of freedom in the usual ANOVA: 

= v& 

®Reg 



Giving cj) a CAR prior is equivalent to giving a N(0,tD) prior; the latter 
are the "constraint cases" in HCSC's SANOVA structure. The precision 
Dnn = for the overall level is equivalent to a flat prior on @gm, though 
@GM could alternatively have a normal prior with mean zero and finite 
variance. If G > 1, the CAR prior also puts zero precision on G — 1 contrasts 
in <f>, which are contrasts in the levels of the G islands [Hodges, Carlin and 
Fan (2003)]. 



2.2. SANOVA as a competitor to MCAR. 



2.2.1. Multivariate conditionally autoregressive (MCAR) models. With 
multiple diseases, we have unknown faj corresponding to region i and dis- 
ease j, where i = 1, . . . , N and j = 1, . . . , n. Letting fibea common precision 
matrix (i.e., inverse of the covariance matrix) representing correlations be- 
tween the diseases in a given region, MCAR distributions arise through 
conditional specifications for (f) i = ((f>n, . . . , (f>i n )'' 

(8) 0il{«Mi'* ~ mvn( — £ — n- 1 ) . 

These conditional distributions yield a joint distribution for cj) = ■ ■ ■ , 0jv)' : 

(9) f{<t>\n) oc \\n\\( N - G V 2 exp {-\4>'{Q ® 0)0), 

where Q is defined as in Section 2.1.2 and again (9) is an improper den- 
sity when a = 1. However, as for the univariate CAR, this yields proper 
posteriors in conjunction with a proper likelihood. The specification above 
is a "separable" dispersion structure, that is, the covariances between the 
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diseases are invariant across regions. This may seem restrictive, but relax- 
ing this restriction gives even more complex dispersion structures [see Jin, 
Banerjee and Carlin (2007) and references therein]. As mentioned earlier, our 
focus is to retain the model's simplicity without compromising the primary 
inferential goals. We propose to do this using SANOVA and will compare it 
with the separable MCAR only. 



2.2.2. SANOVA with Minnesota counties as one factor. We now de- 
scribe the SANOVA model using the Minnesota 3-cancer dataset. Con- 
sider the Minnesota map with N = 87 counties, and suppose each county 
has counts for n = 3 cancers. County i has an n- vector of parameters de- 



scribing the n cancers, (f) i 



, define the Nn vector as 



(f> = (<j>i, (f}' 2 , . . . , (^'n)' ■ For now, we are vague about the specific interpreta- 
tion of 4>ij\ the following description applies to any kind of data. Assume 
the N x N matrix Q describes neighbor pairs among counties as before. 
The SANOVA model for this problem is a 2-way ANOVA with factors can- 
cer ("CA," n levels) and county ("CO," iV levels) and no replication. As 
in Section 2.1.1, we model <j) with a saturated linear model and put the 
grand mean and the main effects in their traditional positions as in ANOVA 
(matrix dimensions and definitions appear below the equation): 



cj ) =[cf>' 1 ,cf ) ' 2 ,...,ct ) / N ]' = [A 1 \A 2 ]e 



(10) 



'Nn 



LiVn 



j^l N ®H c y 



Grand mean 
Nnxl 



Cancer 
main effect 
Nnx{n— 1) 



■CA 



County 
main effect 
Nnx(N-l) 



Cancer X County 
interaction 
Nnx(N-l)(n-l) 



©GM 
&CA 
X &CO ' 
-®COxCA- 

where Hca is an n x (n — 1) matrix whose columns are contrasts among 
cancers, so l' n HcA = 0^-1' ana - H' ca Hca = In-i] Hqa is the jth column 
of Hca j an d is V without its Nth column J=ljy, so it has N — 1 

columns, each a contrast among counties, that is, V N V^ = 0^ r _ 1 , and 
y(-)/y(-) =I N _ lm The column labeled "Grand mean" corresponds to the 
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ANOVA's grand mean and has parameter ©gm! the other blocks of columns 
labeled as main effects and interactions correspond to the analogous AN OVA 
effects and to their respective parameters ®ca,®co,®COxCA- Defining 
prior distributions on completes the SANOVA specification. We put in- 
dependent flat priors (normal with large variance) on &gm and &ca, which 
are, therefore, not smoothed. This is equivalent to putting a flat prior on 
each of the n cancer-specific means. To specify the smoothing priors, de- 
fine Hq A = -7= l n . Let the county main effect parameter &co have prior 

®CO ~ ^Yzv-i(O) r oD^), where Z?(~) corresponds to V^~\ that is, D with- 
out its iVth row and column, To > is unknown and tqD^~' is a precision 
matrix. Similarly, let the jth group of columns in the cancer-by-county in- 



teraction, V<-> ®H$ A% have prior 



1 COxCA 



N N -i(0,TjD(-)), for Tj >0 



(7) 

unknown. Each of the priors on &co and the ®coxCA * s a CAR prior; the 
overall level of each CAR, with prior precision zero, has been included in 
the grand mean and cancer main effects. 

To compare this to the MCAR model, use SANOVA's priors on to 
produce a marginal prior for cf) comparable to the MCAR's prior on cf) (Sec- 
tion 2.2.1); in other words, integrate @co and ®coxCA out of the foregoing 
setup. A priori, 



(11) 



K 







©CO 

(1) 

COxCA 



\ 



V R (n_1) / 
\ CO xCA / 



has precision Q®{H A +S} dia,g{Tj)H A + ^'), where K is the columns of the design 
matrix for the county main effects and cancer-by-county interactions — the 
right-most n(N — 1) columns in equation (10)'s design matrix — and H A + ^ = 
(-j=1_ n \H ca) is an orthogonal matrix. Appendix A in Zhang, Hodges and 
Banerjee (2009) gives a proof. 



2.2.3. Comparing SANOVA vs MCAR. Defining cf) as in Sections 2.2.1 
and 2.2.2, consider the MCAR prior for cf), with within-county precision 
matrix fi. Let have spectral decomposition VciDqV^, where isnxn 
diagonal and Vq is n x n orthogonal. Then the prior precision of cf) is Q <S> 
(VqDqVq), where Q is the known neighbor relations matrix and Vh and Dq 
are unknown. Comparing MCAR to SANOVA, the prior precision matrices 
for the vector cf) are as in Figure 1. SANOVA is clearly a special case of 



MCAR in which H A + ^ is known. Also, as described so far, H K A ' has one 
column proportional to l n with the other columns being contrasts, while 
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MCAR: 



SANOVA: 



Q (VqDqVq) 
t t t 



unknown 



Q ® (H/+'>diag(T j )H/ + y) 



known- 



unknown^ 



"-unknown 



Fig. 1. Comparing prior precision matrices for <fi in MCAR and SANOVA. 



MCAR avoids this restriction. MCAR is thus more flexible, while SANOVA 
is simpler, presumably making it better identified and easier to set priors for. 
MCAR should have its biggest advantage over SANOVA when the "true" 
Vq is not like for any specification of the smoothing precisions tj. 

However, because data sets often have modest information about higher- 
level variances, it may be that using the wrong usually has little effect 
on the analysis. In other words, SANOVA's performance may be relatively 
stable despite having to specify \ while MCAR may be more sensitive 
to O's prior. 



2.3. Setting priors in MCAR and SANOVA. 



2.3.1. Priors in SANOVA. For the case of normal errors, based on equa- 
tions (1) and (10), setting priors for @,Tj, t]q completes a Bayesian specifica- 
tion. Since r and r/o are precision parameters, one possible prior is Gamma; 
this paper uses a Gamma with mean 1 and variance 10. As mentioned, the 
grand mean and cancer main effects 61,62,83 have flat priors with tt(8) oc 1, 
though they could have proper informative priors. The priors for 64, . . . are 
set according to the SANOVA structure as in Section 2.2.2. We ran chains 
drawing in the order 6, r and 770 [Appendix B in Zhang, Hodges and Banerjee 
(2009) gives details]. Hodges, Carlin and Fan (2007) also considered priors on 
the degrees of freedom in the fitted model, some conditioned so the degrees 
of freedom in the model's fit were fixed at a certain degree of smoothness. 
The present paper emphasizes comparing MCAR and SANOVA, so we do 
not consider such priors. For the case of Poisson errors, we use a normal 
prior with mean and variance 10 6 for the grand mean and cancer main 
effects 61,62,63. The other 6{S are given normal CAR priors as discussed in 
Section 2.2.2. For the prior on the smoothing precisions Tj, we use Gamma 
with mean 1 and variance 10. To reduce high posterior correlations among 
the 6s, we used a transformation during MCMC; Appendix C in Zhang, 
Hodges and Banerjee (2009) gives details. 
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2.3.2. Priors in MCAR. MCAR models were fitted in WinBUGS. For 
the normal-error case, we used this model and parameterization: 



i = 1, . . . ,N; j = 1, . . . ,ra, where rjo has a gamma prior with mean 1 and 
variance 10 as for SANOVA. To satisfy WinBUGS 's constraint that ^ Sij = 
0, we add cancer-specific intercepts f3j. We give j3j a flat prior and for S, the 
spatial random effects, we use an intrinsic multivariate CAR prior. Similarly, 
in the Poisson case 



where Eij is an offset. Prior settings for {3j and Sij are as in the nor- 
mal case. For MCAR priors, the within-county precision matrix Q needs 
a prior; a Wishart distribution is an obvious choice. If f2 ~ Wishart(-R, v), 
then E(Q) = vRr 1 . We want a "vague" Wishart prior; usually v = n is used 
but little is known about how to specify R. Thus, we considered three dif- 
ferent Rs, each proportional to the identity matrix. One of these priors sets 
i?'s diagonal entries to Ra = 0.002, close to the setting used in an example 
in the GeoBUGS manual (oral cavity cancer and lung cancer in West York- 
shire). The other two -Rs are the identity matrix and 200 times the identity. 
For the special case n = 1, where the Wishart reduces to a Gamma, these 
Wisharts are T(0.5, 0.001), T(0.5,0.5) and r(0.5,100), respectively. 

3. Simulation experiment. For this simulation experiment, artificial data 
were simulated from the model used in SANOVA with a spatial factor, as 
described in Section 2.2.2. Three different types of Bayesian analysis were 
applied to the simulated data: SANOVA with the same used to gener- 
ate the simulated data (called "SANOVA correct"); SANOVA with incorrect 
H { / ] ; and MCAR. SANOVA correct is a theoretical best possible analysis 
in that it takes as known things that MCAR estimates, that is, it uses ad- 
ditional correct information. SANOVA correct cannot be used in practice, 
of course, because the true is not known. MCAR vs SANOVA with 

incorrect is the comparison relevant to practice, and comparing them 

to SANOVA correct shows how much each method pays for its "deficiency" 
relative to SANOVA correct. Obviously it is not enough to test the SANOVA 
model using only data generated from a similar SANOVA model. To avoid 
needless computing and facilitate comparisons, instead of generating data 
from an MCAR model and fitting a SANOVA model as specified above, we 
use a trick that is equivalent to this. Section 3.1.2 gives the details. 



(12) 





(13) 



Yij ~ Poisson 
log(mj) = log(Eij) + /3j + S^, 
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3.1. Design of the simulation experiment. We simulated both normally- 
distributed and Poisson-distributed data. For both types of data, we consid- 
ered two different true sets of smoothing parameters r = r/r/o or r (Table 
1). For the normal data, we considered t/tjo, since this ratio determines 
smoothing in normal models, and we also considered two error precisions rjo 
(Table 1). 

3.1.1. Generating the simulated data sets. To generate data from the 
SANOVA model, we need to define the true H^~\ Let 



HA, 





V 




We used HA\ as the correct its columns are scaled to have length 1. 

Given and with known, one draw of and e produces a draw of 
Xr>® + e, therefore a draw of y. In the simulation, we let the grand mean 
and main effects, which are not smoothed, have true value 5. Each observa- 
tion is simulated from a 3 x 20 factorial design, where 3 is the number of 
cancers and 20 is the number of counties. We used the 20 counties in the 
right lower corner of Minnesota's map, with their actual neighbor relations. 
Thus, the dimension of each artificial data set is 60. The simulation experi- 
ment is a repeated-measures design, in which a "subject" s in the design is 
a draw of (<$^,"y^), referring to equation (3), where 8^}_ s = 5 and S^Iqq ~ 
^57(0,13 <g> D(~>) specify and ~ A^6o(0)^6o) gives e. For the normal- 
errors case, 100 such "subjects" were generated. Given a design cell in the 
simulation experiment with r = (a, b, c) and t]q = d, the artificial data set for 
subject a is yW=X D diag(l' 3 , -*-l' 19 , -i.l' 19> -i.l' 19 )tfW + -^ 7 W. All fac- 
tors of the simulation experiment were applied to each of the 100 "subjects." 
For the normally-distributed data, the simulation experiment had these fac- 
tors: (a) the true (t /t?o, n/'70,^/ 7 ») : ( 100 ' 100 ' ai ) or C - 1 ' 100 ' 0A )' ( b ) the 

Table 1 

Experimental conditions in the simulation experiments 
Error distribution r} (to/j7 ,ti/t7 ,T2/t7 )|(to,ti,T2) Data name 



Normal 


1 


(100, 100, 0.1) 


Datal 




1 


(0.1, 100, 0.1) 


Data2 




10 


(100, 100, 0.1) 


Data3 




10 


(0.1, 100, 0.1) 


Data4 


Poisson 


NA 


(100, 100, 0.1) 


Data5 




NA 


(0.1, 100, 0.1) 


Data6 
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true error precision t]q: 1 or 10; and (c) six statistical methods, described 
below in Section 3.1.2. Each design cell described in Table 1 thus had 100 
simulated data sets. Similarly, for the Poisson-data experiment, another 100 
"subjects" were generated, but now there is no 7^). Thus, each "subject" 
s is a vector S^ s \ where 5^ is as described above. For the design cell with 
r = (100,100,0.1), the artificial data for subject s is ~ Poisson(^( s )), 
where log(/xW) = \og(E)+X D diag(l' 3 , ^1' 19 , ^1' 19 , -^=l' 19 )tfW. In the sim- 
ulation experiment, we use "internal standardization" of the Minnesota 3- 
cancer data to supply the expected numbers of cancers Eij. Among the 20 
extracted counties, Hennepin county has the largest average population over 
11 years, about 1.1 million; its cancer counts are 5294, 119 and 439 for lung, 
larynx and esophagus respectively. Faribault county has the smallest average 
population, 16,501, with cancer counts 110, 7 and 13 respectively. The 
have ranges 80 to 5275, 2 to 113 and 7 to 449 for lung, larynx and esoph- 
agus cancer respectively. For the Poisson data, the simulation experiment 
had these factors: (a) the true tq,t\,T2' (100,100,0.1) or (0.1,100,0.1); and 
(b) six statistical methods described below in Section 3.1.2. Again, each of 
the two design cells in Table 1 had 100 simulated data sets. 



3.1.2. The six methods (procedures). For each simulated data set, we 
did a Bayesian analysis for each of six different models described in Table 2. 
The six models are: SANOVA with the correct HA 1 ; SANOVA with a 

somewhat incorrect , HA 2 given below; a variant SANOVA with a very 
incorrect H A + \ H AM given below; MCAR with Ru = 0.002; MCAR with 
Ru = 1; and MCAR with Ru = 200 (see Section 2.3.2). HA 2 and H AM are 

/ll l\/73 ° 
HA 2 =[ 1 -2 ^ 

V 1 -i/\o 




Table 2 

The six statistical methods used in the simulation experiment 



Procedure 



SANOVA with correct H { /> 

SANOVA with incorrect H A ^ 

Variant SANOVA with Ham 

MCAR 

MCAR 

MCAR 



Prior 



rjd,Ti~r(0.1,0.1) for 3 = 0,1, 2 
ij0,Tj~r(0.1,0.lj for 3= 0,1, 2 
ijo,r i ~r(0.1,0.1) fori = 0,1, 2 
n ~ Wishart(i?, 3), R = 0.002J 3 

n ~ Wishart(7?, 3), R = h 
n ~ Wishart(7i, 3), R= 200/ 3 
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/ 0.56 -0.64 -0.52\ 
H AM = -0.53 -0.77 0.36 . 
\-0.63 0.07 -0.77/ 

The incorrect HA2 has the same first column (grand mean) as the correct 
HA\, so it differs from the correct HAi, though less than it might. As noted 
above, we need to see how the SANOVA model performs for data generated 
from an MCAR model in which Vh from Figure 1 does not have a column 
proportional to l n . To do this without needless computing, we used a trick: 
we used the data generated from a SANOVA model with HA\ and fit the 
variant SANOVA mentioned above, in which H^' is replaced by the or- 
thogonal matrix Ham with no column proportional to l n , chosen to be very 
different from HA\. For normal errors (Datal through Data4), this is pre- 
cisely equivalent to fitting a SANOVA with Ha — HA\ to data generated 
from an MCAR model with Vn = BHA 1 , for B = HAiH^ M , that is, 

0.43 -0.74 -0.52" 
-0.13 -0.63 0.77 
-0.89 -0.26 -0.37 

(to 2 decimal places). For Poisson errors (Data5, Data6), the equivalence 
is no longer precise but the divergence of fitted SANOVA [using H A + ^ = 

Ham] and generated data [using H^ = HA±] is quite similar. Finally, we 
considered three priors for MCAR because little is known about how to set 
this prior and we did not want to hobble MCAR with an ill-chosen prior. 
For the SANOVA and variant SANOVA analyses, we gave tj a T(0. 1,0.1) 
prior with mean 1 and variance 10 for both the normal data and the Poisson 
data. 



(14) 



3.2. Outcome measures. To compare the six different methods for nor- 
mal and Poisson data, we consider three criteria. The first is average mean 
squared error (AMSE). For each of the 60 (X^Q)-, the mean squared error 
is defined as the average squared error over the 100 simulated data sets. 
AMSE for each design cell in the simulation experiment is defined as the 
average of mean squared error over the 60 {X^®)^. Thus, for the design 
cell labeled DataK in Table 1, define 

L N n 

(15) AMSE K = r EEEK^ S )S " (XDS^/Nn, 

d=i i=i j=i 

where L = 100, N = 20, n = 3, if = 1, . . . , 4 for Normal, K = 5, 6 for Pois- 
son, is the true value and is the posterior median of 0. For each 
design cell (K), the Monte Carlo standard error for AMSE is (100) -0 ' 5 
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Table 3 

AMSE for simulated normal and Poisson data 



Normal-error model Poission-error model 



Model 


Datal 


Data2 


Data3 


Data4 


Data5 


Data6 


SANOVA with HAi 


0.34 


0.60 


0.04 


0.06 


0.02 


0.04 


SANOVA with HA 2 


0.47 


0.84 


0.05 


0.07 


0.02 


0.14 


SANOVA with HAam 


0.48 


0.74 


0.05 


0.06 


0.03 


0.11 


MCAR with R u = 0.002 


0.66 


1.88 


0.04 


0.13 


0.02 


0.04 


MCAR with Ru = 1 


0.36 


0.84 


0.04 


0.06 


0.02 


0.06 


MCAR with Ru = 200 


0.93 


0.92 


0.09 


0.09 


0.24 


0.36 



times the standard deviation, across DataK's 100 simulated data sets, of 
EiIiEj=iP^©)J - {X D ®) d i3 ] 2 /Nn. The second criterion is the bias of 
Xe>@. For each of DataK's 100 simulated data sets, first compute posterior 
medians of (XdO) 1 1 , . . . , (Xd®) 20 3 , then average each of those posterior 
medians across the 100 simulated data sets. From this average, subtract 
the true (Xf)@) i jS to give the estimated bias for each of the 60 (Xo®) i jS. 
MBIAS is defined as the 2.5th, 50th and 97.5th percentiles of the 60 esti- 
mated biases. More explicitly, for design cell DataK, MBIAS is 

MBIAS K = 2.5th, 50th, 97.5th percentiles of 

(16) 

(iE(^© d -^))- 

Finally, the coverage rate of Bayesian 95% equal-tailed posterior intervals, 
"PI rate," is the average coverage rate for the 60 individual (X^Q)^. 

3.3. Markov chain Monte Carlo specifics. While the MCAR models were 
implemented in WinBUGS, our SANOVA implementations were coded in R 
and run on Unix. The different architectures do not permit a fair compar- 
ison between the run times of SANOVA and MCAR. However, the SANOVA 
models have lower computational complexity than the MCAR models: MCAR 
demands a spectral decomposition in every iteration, while SANOVA does 
not. For each of our models, we ran three parallel MCMC chains for 10,000 
iterations. The CODA package in R (www.r-project.org) was used to di- 
agnose convergence by monitoring mixing using Gelman-Rubin diagnostics 
and autocorrelations [e.g., Gelman et al. (2004), Section 11.6]. Sufficient 
mixing was seen within 500 iterations for the SANOVA models, while 200 
iterations typically revealed the same for the MCAR models; we retained 
8000 x 3 samples for the posterior analysis. 
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3.4. Results. Table 3 and Figures 2 and 3 show the simulation experi- 
ment's results. Table 3 shows AMSE; for all methods and design cells, the 
standard Monte Carlo errors of AMSE are small, less than 0.07, 0.005 and 
0.025 for Datal/Data2, Data3/Data4 and Data5/Data6 respectively. Figure 
2 shows MBIAS, where the middle symbols represent the median bias and 
the line segments represent the 2.5th and 97.5th percentiles. Figure 3 shows 
coverage of the 95% posterior intervals. Denote SANOVA with the correct 
H { + ] {HA]) as "SANOVA correct," SANOVA with HA 2 as "SANOVA in- 
correct," the variant SANOVA with H AM as "SANOVA variant," MCAR 
with Ru = 0.002 as "MCAR0.002" and so on. 

3.4.1. As expected, SANOVA with correct H^ performs best. For nor- 
mal data, SANOVA correct has the smallest AMSE for all true r/o and t 
(Table 3). The advantage is larger in Datal and Data2 where the error pre- 
cision r/o is 1 than in Data3 and Data4 where tjq is 10 (i.e., error variation is 
smaller). For Poisson data, SANOVA correct also has the smallest AMSE. 
Considering MBIAS (Figure 2), SANOVA correct has the narrowest MBIAS 
intervals for all cases. In Figure 3, the posterior coverage for SANOVA cor- 
rect is nearly nominal. As expected, then, SANOVA correct performs best 
among the six methods. 

3.4.2. SANOVA with incorrect HA 2 and Ham perform very well. Table 3 
shows that, for normal data, both SANOVA incorrect and SANOVA variant 
have smaller AMSEs than MCAR200 and MCAR0.002, and AMSEs at worst 
close to MCARl's. For Poisson data, Table 3 shows that MCAR0.002 and 
MCAR1 do somewhat better than SANOVA incorrect and variant SANOVA. 



-e- SANOVA with correct HA1 
a- SANOVA with incorrect HA2 
-V- variant SANOVA with H_AM 
-+- MCAR with Rii=0.002 
-x- MCAR with Rii=1 
MCAR with Rii=200 



-G- 



(a) Normal (b) Poisson 

Fig. 2. MBIAS for simulated normal and Poisson data. 
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SANOVA with correct HA1 
SANOVA with incorrect HA2 
O variant SANOVA with H AM 
-O- MCAR with Rii=0.002 
v- MCAR with Rik1 
-0- MCAR with Rii=200 



(a) Normal (b) Poisson 

Fig. 3. PI rate for simulated normal and Poisson data. 



Considering MBIAS in normal data [Figure 2(a)], the width of the 95% 
MBIAS intervals for SANOVA incorrect are the same as or smaller than 
for all three MCAR procedures. Similarly, SANOVA variant has MBIAS 
intervals better than MCAR0.002 and MCAR200 and almost as good as 
MCAR1. Figure 2(b) for Poisson data shows SANOVA correct, MCAR0.002 
and MCAR1 have similar MBIAS intervals. SANOVA variant in Data5 and 
SANOVA incorrect in Data6 show the worst performance for MBIAS apart 
from MCAR200, whose MBIAS interval is much the widest. Figure 3(a) 
shows that, for normal data, interval coverage for SANOVA incorrect and 
SANOVA variant is very close to nominal. It appears that the specific value 
of has little effect on PI coverage rate for the cases considered here. 

Apart from MCAR200 for Datal/Data2 and MCAR0.002 for Datal through 
Data4, which show low coverage, all the other methods have coverage rates 
greater than 90% for normal data, most close to 95%. For Data3 and Data4, 
PI rates for MCAR200 reach above 99%. For Poisson data, the PI rates for 
SANOVA incorrect and SANOVA variant are close to nominal and better 
than MCAR0.002 and MCAR200. In particular, all SANOVAs have the 
closest to nominal coverage rates for both normal and Poisson data, which 
again shows the stability of SANOVA under different settings. 

3.4.3. MCAR is sensitive to the prior on £1. To fairly compare SANOVA 
and MCAR, we considered MCAR under three different prior settings. For 
normal data, MCAR1 has the smallest AMSEs and narrowest MBIAS in- 
tervals among the MCARs considered, while MCAR0.002 has the largest 
and widest, respectively. For Poisson data, however, MCAR0.002 has the 
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best AMSE and MBIAS among the MCARs. MCAR200 performs poorly 
for both Normal and Poisson. The coverage rates in Figure 3 show similar 
comparisons. These results imply that the prior matters for MCAR: no sin- 
gle prior was always best. By comparison, SANOVA seems more robust, at 
least for the cases considered. 

3.5. Summary. As expected, SANOVA correct had the best performance 
because it uses more correct information. For normal data, SANOVA incor- 
rect and SANOVA variant had similar AMSEs, better than two of the three 
MCARs for the data sets considered. For Poisson data, SANOVA incorrect 
and SANOVA variant had AMSEs as good as those of MCAR0.002 and 
MCAR1 for Data5 and somewhat worse for Data6, while showing nearly 
nominal coverage rates in all cases and less tendency to bias than MCAR 
in most cases. Replacing the r(0. 1,0.1) prior for r with r(0. 001, 0.001) left 
AMSE and MBIAS almost unchanged and coverage rates a bit worse (data 
not shown). MCAR, on the other hand, seems more sensitive to the prior 
on J7. MCAR0.002 tends to smooth more than MCAR1, more so in normal 
models where the prior is more influential than in Poisson models. (The lat- 
ter is true because data give more information about means than variances, 
and the Poisson model's error variance is the same as its mean, while the nor- 
mal model's is not.) For the normal data, MCAR0.002's tendency to extra 
shrinkage appears to make it oversmooth and perform poorly for Data2 and 
Data4, where the truth is least smooth. For the Poisson data, MCAR0.002 
and MCAR1 give results similar to each other and somewhat better than 
the SANOVAs except for interval coverage. Therefore, SANOVA, with stable 
results under different H^' and with parameters that are easier to under- 
stand and interpret, may be a good competitor to MCAR in multivariate 
spatial smoothing. 

4. Example: Minnesota 3-cancer data. Researchers in different fields have 
illustrated that accounting for spatial correlation could provide insights that 
would have been overlooked otherwise, while failure to account for spatial 
association could potentially lead to spurious and sometimes misleading re- 
sults [see, e.g., Turechek and Madden (2002), Ramsay, Burnett and Krewski 
(2003), Lichstein et al. (2002)]. Among the widely investigated diseases are 
the different types of cancers. We applied SANOVA and MCAR to a cancer- 
surveillance data set describing total incidence counts of 3 cancers (lung, 
larynx, esophagus) in Minnesota's 87 counties for the years 1990 to 2000 
inclusive. Minnesota's geography and history make it plausible that disease 
incidence would show spatial association. Three major North American land 
forms meet in Minnesota: the Canadian Shield to the north, the Great Plains 
to the west, and the eastern mixed forest to the southeast. Each of these re- 
gions is distinctive in both its terrain and its predominant economic activity: 
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mining and outdoors tourism in the mountainous north, highly mechanized 
crop cultivation in the west, and dairy farming in the southeast. The differ- 
ent regions were also settled by somewhat different groups of in-migrants, 
for example, disproportionately many Scandinavians in the north. These 
factors imply spatial association in occupational hazards as well as culture, 
weather, and access to health care especially in the thinly-populated north, 
which might be expected to produce spatial association in diseases. With 
multiple cancers one obvious option is to fit a separate univariate model for 
each cancer. But diseases may share the same spatially distributed risk fac- 
tors, or the presence of one disease might encourage or inhibit the presence 
of another in a region, for example, larynx and esophagus cancer have been 
shown to be closely related spatially [Baron et al. (1993)]. Thus, we may 
need to account for dependence among the different cancers while main- 
taining spatial dependence between sites. Although the data set has counts 
broken out by age groups, for the present purpose we ignore age standard- 
ization and just consider total counts for each cancer. Age standardization 
would affect only the expected cancer counts Eij, while other covariates 
could be added to either SANOVA or MCAR as unsmoothed fixed effects 
(i.e., in the A\ design matrix). Given the population and disease count of 
each county, we estimated the expected disease count for each cancer in each 
county using the Poisson model. Denote the 87 x 3 counts as yi 5 i, . . . ,2/87,3! 
then the model is 

V%j\Vnj ~ Poisson(fe), 

(17) 

log(/iy) = log(-Eij) + (X D Q)ij, 
where Xe>& is the SANOVA structure and & has priors as in Section 2.2.2. 

y . o 

For disease j in county i, = Pi yi p , where Oij is the disease count for 
county i and disease j and Pi is county i's population. For the SANOVA de- 
sign matrix, we consider HA\ and HA2 from the simulation experiment, 
though now neither is known to be correct. We also consider a variant 
SANOVA analysis using estimated from the MCAR1 model, to test 

the stability of the SANOVA results. Appendix D in Zhang, Hodges and 
Banerjee (2009) describes the latter analysis. Figures 4 to 6 show the data 
and results for MCAR1 and SANOVA with HA\. In each figure, the upper 
left plot shows the observed y^/Ey; the two lower plots show the posterior 
median of fiij/Eij for MCAR1 and SANOVA with HA\. Lung cancer counts 
tended to be high and thus were not smoothed much by any method, while 
counts of the other cancers were much lower and thus smoothed consider- 
ably more (see also Figure 7). Since SANOVA with HAi, HA2 and estimated 
gave very similar results, only those for HA\ are shown. Results for 
MCAR0.002 are similar to those for MCAR1, so they are omitted. As ex- 
pected, MCAR200 shows the least shrinkage among the three MCARs and 
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gives some odd Hij/Eij. To compare models, we calculated the Deviance In- 
formation Criterion [DIC; Spiegelhalter et al. (2002)]. To define DIC, define 
the deviance D(9) = — 21og/(y|0) + 21og/i(y), where 6 is the parameter 
vector in the likelihood and h(y) is a function of the data. Since h does not 
affect model comparison, we set log h(y) to 0. Let 6 be the posterior mean of 
and D the posterior expectation of D{6). Then define pr> = D(9) — D(0) 
to be a measure of model complexity and define DIC = D + pp. Table 4 
shows D,po and DIC for nine analyses, SANOVA with 3 different 
MCAR with 3 different priors for O, and 3 fits of univariate CAR mod- 
els to the individual diseases, discussed below. Considering D, the three 
SANOVAs and MCAR1 are similar; Figures 4 to 6 show the fits are indeed 
similar. Figure 7 reinforces this point, showing that MCAR1 and SANOVA 
with HA\ induce similar smoothing for the three cancers. SANOVA with 
estimated from MCAR has the smallest D (1458), though its model 
complexity penalty (p D = 103) is higher than MCAR0.002's (p D = 79). De- 
spite having the second worst fit (D), MCAR0.002 has the best DIC, and 
the three SANOVAs have DICs much closer to MCAR0.002's than to the 
other MCARs. Generally, all SANOVA models have similar D (s» 1460) and 
DIC (~ 1562), while MCAR results are sensitive to fi's prior, consistent 
with the simulation experiment. For comparison, we fit separate univariate 
CAR models to the three diseases considering three different priors for the 
smoothing precision, r~Gamma(a,o) for a = 0.001, 1 and 1000. For each 
prior, we added up D, pu and DIC for three diseases (see Table 4). With 
a = 0.001 and 1, we obtained D's (1461 and 1453 respectively) competi- 
tive with SANOVA, MCAR0.002 and MCAR1 but with considerably greater 
complexity penalties (141 and 149 respectively) and thus DICs slightly larger 
than 1600. For a = 1000, we obtained an even lower D (1432), but an in- 
creased penalty (180) resulted in a poorer DIC score. Figure 7 shows fitted 
values for CAR1, which were smoothed like MCAR1 and SANOVA for lung 
and esophagus cancers but smoothed rather more for larynx cancer. Over- 
all, these results reflect some gain in performance from accounting for the 
space-cancer interactions /associations. 

To further examine the smoothing under SANOVA, Figure 8 shows sep- 
arate maps for the county main effect and interactions from the SANOVA 
fit with HA\. The upper left plot is the cancer main effect, the mean of 
the three cancers; the lower left plot is the comparison of lung versus av- 
erage of larynx and esophagus; the lower right plot is the comparison of 
larynx versus esophagus. All values are on the same scale as yij/Eij in Fig- 
ures 4 to 6 and use the same legend. The two interaction contrasts are 
smoothed much more than the county main effect, agreeing with previ- 
ous research that larynx and esophagus cancer are closely related spatially 
[Baron et al. (1993)]. To see whether the interactions are necessary, we 
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(b) fitted: SANOVA with HAi (c) fitted: MCAR with R u =1 



Fig. 4. Lung cancer data and fitted values. 

fit a SANOVA model (using HA\) without the interactions. As expected, 
model complexity decreased (pd = 77), while D increased slightly, so DIC 
became 1558, a bit better than SANOVA with interactions. Now consider 
the posterior of the MCAR's precision matrix £1. The posterior mean of 
Q is much larger for MCAR0.002 than MCAR200; the diagonal elements 
are larger by 4 to 5 orders of magnitude. This may explain the poor cov- 
erage for MCAR0.002 in the simulation. Further, consider the correlation 
matrix arising from the inverse of Si's posterior mean. As the diagonals 
of R change from 0.002 to 200, the correlation between any two cancers 
decreases and the complexity penalty po increases. By comparison, the 
three SANOVAs have similar model fits and complexity penalties, lead- 
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(a) data: yij/Eij 
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(b) fitted: SANOVA with HA X (c) fitted: MCAR. with R„ =1 

Fig. 5. Larynx cancer data and fitted values. 

ing to similar DICs. So again, in this sense SANOVA shows greater sta- 
bility. 

5. Discussion and future work. We used SANOVA to do spatial smooth- 
ing and compared it with the much more complex MCAR model. For the 
cases considered here, we found SANOVA with spatial smoothing to be an 
excellent competitor to MCAR. It yielded essentially indistinguishable in- 
ference, while being easier to fit and interpret. In the SANOVA model, 

is assumed known. For most of the SANOVA models considered, s 
first column was fixed to represent the average over diseases, while other 
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(b) fitted: SANOVA with HAi (c) fitted: MCAR with R H =1 



Fig. 6. Esophagus cancer data and fitted values. 

columns were orthogonal to the first column. Alternatively, could be 

treated as unknown and estimated as part of the analysis. With this ex- 
tension, SANOVA with spatial effects is a reparameterization of the MCAR 
model and gains the MCAR model's flexibility at the price of increased 
complexity. This extension would be nontrivial, involving sampling from the 
space of orthogonal matrices while avoiding identification problems arising 
from, for example, permuting columns of Other covariates can be 

added to a spatial SANOVA. Although (10) is a saturated model, spatial 
smoothing "leaves room" for other covariates. Such models would suffer from 
collinearity of the CAR random effects and the fixed effects, as discussed by 
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Fig. 7. Comparing data and fitted values for each cancer. The "Data" panel shows the 
density for Uij/Eij, while the other three panels show the posterior median of inj/Eij for 
univariate CAR, SANOVA and MCAR. 



Reich, Hodges and Zadnik (2006), who gave a variant analysis that avoids 
the collinearity. For data sets with spatial and temporal aspects, for exam- 
ple, the 11 years in the Minnesota 3-cancer data, interest may lie in the 
counts' spatial pattern and in their changes over time. By adding a time 
effect, SANOVA can be extended to a spatiotemporal model. Besides spa- 
tial and temporal main effects, their interactions can also be included and 
smoothed. There are many modeling choices; the simplest model is an ad- 
ditive model without space-time interactions, where the spatial effect has a 
CAR model and the time effect a random walk, which is a simple CAR. But 
many other choices are possible. We have examined intrinsic CAR models, 
where Qij = — 1 if region i and region j are connected. SANOVA with spatial 
smoothing could be extended to more general CAR models. Banerjee, Carlin 
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Table 4 
Model comparison using DIC 



Model 


D 


Pd 


DIC 


SANOVA with HAi 


1461 


103 


1564 


SANOVA with HA 2 


1463 


102 


1565 


SANOVA with H A estimated from MCARl 


1458 


103 


1561 


MCAR0.002 


1476 


79 


1555 


MCARl 


1459 


132 


1591 


MCAR200 


1559 


356 


1915 


CARO.OOl 


1461 


141 


1602 


CARl 


1453 


149 


1602 


CARIOOO 


1432 


180 


1612 



and Gelfand (2004) replaced Q with the matrix D w — pW, where D w is diag- 
onal with the same diagonal as Q and Wij = 1 if region i is connected with 
region j, otherwise Wij = 0. Setting p = 1 gives the intrinsic CAR model 
considered in this paper. For known p, the SANOVA model described here 
is easily extended by replacing Q in Section 2 with D w — pW. However, for 
unknown p, our method cannot be adjusted so easily, because updating p in 
the MCMC would force V and the design matrix to be updated as well, but 
this would change the definition of the parameter 0. Therefore, a different 
approach is needed for unknown p. A different extension of SANOVA would 
be to survival models for areal spatial data [e.g., Li and Ryan (2002), Baner- 
jee, Wall and Carlin (2003), Diva, Dey and Banerjee (2008)]. If the regions 
are considered strata, then random effects corresponding to nearby regions 
might be similar. In other words, we can embed the SANOVA structure in a 
spatial frailty model. For example, the Cox model with SANOVA structure 
for subject j in stratum i is 

(18) h(tij : X i:j ) = h (tij)exp(XijP), 

where X is the design matrix, which may include a spatial effect, a tem- 
poral effect, their interactions and other covariates. Banerjee, Carlin and 
Gelfand (2004) noted that in the CAR model, considering both spatial and 
nonspatial frailties, the frailties are identified only because of the prior, so 
the choice of priors for precisions is very important. Besides the above ex- 
tensions, HCSC introduced tools for normal SANOVA models that can be 
extended to nonnormal SANOVA models. For example, HCSC defined the 
degrees of freedom in a fitted model as a function of the smoothing preci- 
sions. This can be used as a measure of the fit's complexity, or a prior can 
be placed on the degrees of freedom as a way of inducing a prior on the 
unknowns in the variance structure. The latter is under development and 
will be presented soon. 
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(b) (c) 

Fig. 8. SANOVA with HA\: (a.) county main effect; (h) cancer x county interaction 1 
for larynx; (c) cancer x county interaction 2 for esophagus. 
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SUPPLEMENTARY MATERIAL 

Appendices, data and code (DOI: 10.1214/09-AOAS267SUPP; .zip). Our 
supplementary material includes four sections as appendices. In Appendix 
A we present a derivation of the precision matrix of (11). Details of our 
MCMC algorithms can be found in Appendix B. Appendix C discusses the 
mean transformation for the Poisson case, while Appendix D discusses the 



26 



Y. ZHANG, J. S. HODGES AND S. BANERJEE 



estimation of the H\ from MCAR1. In addition, we provide a compressed 
folder containing the data set for our 3-cancer Minnesota example as well 
as an R code example to implement the SANOVA models. 
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